unit etools;

interface
uses common;

procedure node_pair_energy(n1,n2: integer;Var eplus:myreal);
procedure atom_node_energy(po,t : integer;r:v3;Var eplus:myreal);
procedure atom_pair_energy(identic :boolean; r1,r2 :v3; t1,t2 : integer;Var eplus:myreal);
procedure UPPN(Var epot: myreal);
procedure chempot(Var ech :  myreal);

implementation
uses tools,ppot;

procedure node_pair_energy(n1,n2: integer;Var eplus:myreal);
Var
   i,j	       : integer;
   o,p1,p2,p3  : integer;
   jstart      : integer;
   p	       : array[1..3] of integer;
   v,w	       : v3;
   ar,u,ud,udd : myreal;
   IdenticNode : boolean;
   IdenticAtom : boolean;
begin
   eplus:=0.0;
   IdenticNode:=(n1=n2);
   for i:=1 to nan[n1] do if (n_a[n1,i,1]>0) then begin
      {prevent double-counting when n1=n2 ("self" energy)}
      if not(IdenticNode) then jstart:=1 else jstart:=i;
      for j:=jstart to nan[n2] do if (n_a[n2,j,1]>0) then
	 for p1:=-iext[1] to iext[1] do
	    for p2:=-iext[2] to iext[2] do
	       for p3:=-iext[3] to iext[3] do begin
		  p[1]:=p1; p[2]:=p2; p[3]:=p3;
		  IdenticAtom:=((IdenticNode)
				and (i=j) and (p1=0) and (p2=0) and (p3=0));
		  for o:=1 to 3 do v[o]:=bound(n_r[n2,j,o]+rt[n2,o]-n_r[n1,i,o]-rt[n1,o])+p[o];
		  ortho3(v,w);
		  ar:=sqrt(w[1]*w[1]+w[2]*w[2]+w[3]*w[3]);
		  if ((ar<rcut) and not(IdenticAtom)) then begin
		     upp(n_a[n1,i,1],n_a[n2,j,1],ar,u,ud,udd);
		     {writeln(n1:5,n2:5,i:3,j:3,ar:10:3,u:12:6);}
		     eplus:=eplus+u;
		     if IdenticNode and (i=j) then eplus:=eplus-u/2.0;
		  end;
	       end;
   end;
end; { node_pair_energy }
procedure atom_node_energy(po,t	: integer;r:v3;Var eplus:myreal);
Var
   o,p1,p2,p3,j	: integer;
   p		: array[1..3] of integer;
   v,w		: v3;
   ar,u,ud,udd	: myreal;
begin
   eplus:=0.0;
   if t>0 then begin
      for j:=1 to nan[po] do if (n_a[po,j,1]>0) then
	 for p1:=-iext[1] to iext[1] do
	    for p2:=-iext[2] to iext[2] do
	       for p3:=-iext[3] to iext[3] do begin
		  p[1]:=p1; p[2]:=p2; p[3]:=p3;
		  for o:=1 to 3 do
		     v[o]:=bound(n_r[po,j,o]+rt[po,o]-r[o])+p[o];
		  ortho3(v,w);
		  ar:=sqrt(w[1]*w[1]+w[2]*w[2]+w[3]*w[3]);
		  if ar<rcut then begin
		     upp(t,n_a[po,j,1],ar,u,ud,udd);
		     eplus:=eplus+u;
		  end;
	       end;
   end;
end; { atom_node_energy }
procedure atom_pair_energy(identic :boolean; r1,r2 :v3; t1,t2 : integer;Var eplus:myreal);
Var
   o,p1,p2,p3  : integer;
   p	       : array[1..3] of integer;
   v,w	       : v3;
   ar,u,ud,udd : myreal;
   self,ext0   : boolean;
begin
   eplus:=0.0;
   if ((t1>0) and (t2>0)) then begin
      for p1:=-iext[1] to iext[1] do
	 for p2:=-iext[2] to iext[2] do
	    for p3:=-iext[3] to iext[3] do begin
	       ext0:=((p1=0) and (p2=0) and (p3=0));
	       self:=(ext0 and identic);
	       if not(self) then begin
		  p[1]:=p1; p[2]:=p2; p[3]:=p3;
		  for o:=1 to 3 do
		     v[o]:=bound(r2[o]-r1[o])+p[o];
		  ortho3(v,w);
		  ar:=sqrt(w[1]*w[1]+w[2]*w[2]+w[3]*w[3]);
		  if (ar<rcut) then begin
		     upp(t1,t2,ar,u,ud,udd);
		     eplus:=eplus+u;
		  end;
	       end;
	    end;
   end;
end; { atom_pair_energy }
procedure UPPN(Var epot: myreal);
var i,j	: integer;
   u	: myreal;
begin
   epot:=0.0;
   for i:=1 to tnodes do begin
      for j:=i to tnodes do begin
	 node_pair_energy(i,j,u);
	 epot:=epot+u;
      end;
   end;
end; { UPPN }

procedure chempot(Var ech :  myreal);
Var
   i,j : integer;
begin
   ech:=0;
   for i:=1 to tnodes do
      for j:=1 to nan[i] do if (n_a[i,j,1]>0) then begin
	 ech:=ech+chp[n_a[i,j,1]];
      end;
end; { chempot }

end.			  